Departamento de Física - Faculdade de Ciências e Tecnologia da Universidade de Coimbra

Física Computacional - Ficha 4 - Sistemas de equações Lineares

Rafael Isaque Santos - 2012144694 - Licenciatura em Física

1 - Resolução de um sistema de equações lineares $Ax = b$ pelo método de eliminação de Gauss

In [1]:
import numpy as np
Método de Eliminação de Gauss (sem pivoting)

In [2]:
def gauss_elim_nopiv(n, matrix, b):
    origin_m = np.copy(matrix)
    origin_b = np.copy(b)
    aux_m = np.array([[0.]*n]*n)
    for k in range(0, n):
        # print(origin_m, origin_b)
        for i in range(k+1, n):
            aux_m[i, k] = -origin_m[i, k] / origin_m[k, k]
            origin_m[i, k] = 0
            origin_b[i] = origin_b[i] + aux_m[i, k] * origin_b[k]
            for j in range(k + 1, n):
                origin_m[i, j] = origin_m[i, j] + aux_m[i, k] * origin_m[k, j]

    determ = origin_m[0, 0]
    for d in range(1, n):
        determ *= origin_m[d, d]
    return origin_m, aux_m, origin_b, determ
Backward Substitution

In [3]:
def backward_sub(n, U, c):
    x = np.zeros_like(c)
    x[n-1] = c[n-1]/U[n-1, n-1]
    for k in range(n-2, -1, -1):
        s = 0
        for j in range(k+1, n):
            # print(k, j)
            s += U[k, j] * x[j]
            # print(k, j, U[k, j], x[j], s)
        x[k] = (1/U[k, k]) * (c[k] - s)
        # print(x[k], x)    
    return x
Forward Substitution

In [4]:
def forward_sub(n, L, b):
    y = np.zeros_like(b)
    y[0] = b[1] / L[1, 1]
    for k in range(1, n):
        soma = 0
        for j in range(0, k - 1): soma += L[k, j] * y[j]
        y[k] = (b[k] - soma) / L[k, k]
    return y
Método de Eliminação de Gauss (com pivoting parcial)

In [5]:
def gauss_elim_partial(n, matrix, b):
    origin_m = np.copy(matrix)
    origin_b = np.copy(b)
    
    for k in range(0, n - 1):
        maxi = 0
        km = k
        for k1 in range(k, n):
            if abs(origin_m[k1, k]) > maxi:
                maxi = abs(origin_m[k1, k])
                km = k1
                if km != k:
                    origin_m[k, k:n-1], origin_m[km, k:n-1] = origin_m[km, k:n-1], origin_m[k, k:n-1]
                    origin_b[k], origin_b[km] = origin_b[km], origin_b[k]
    U_A, aux_A, c_A, det_A = gauss_elim_nopiv(n, A, b)
    return U_A, aux_A, c_A, det_A
Funções para resolver sistemas utilizando em primeiro lugar o método sem pivoting e depois com pivoting parcial

In [6]:
def solve_nopiv(matrix, b):
    n = len(matrix)
    u, aux, c, det = gauss_elim_nopiv(n, matrix, b)
    x = backward_sub(n, u, c)
    return x

def solve_partpiv(matrix, b):
    n = len(matrix)
    u, aux, c, det = gauss_elim_partial(n, matrix, b)
    x = backward_sub(n, u, c)
    return x

In [7]:
A = np.array([[1, 2, 1], [2, 3, 3], [-1, -3, 1]])
b = np.array([1, 3, 2])
A, b


Out[7]:
(array([[ 1,  2,  1],
        [ 2,  3,  3],
        [-1, -3,  1]]), array([1, 3, 2]))

In [8]:
gauss_elim_partial(3, A, b)


Out[8]:
(array([[ 1,  2,  1],
        [ 0, -1,  1],
        [ 0,  0,  1]]), array([[ 0.,  0.,  0.],
        [-2.,  0.,  0.],
        [ 1., -1.,  0.]]), array([1, 1, 2]), -1)

In [9]:
gauss_elim_nopiv(3, A, b)


Out[9]:
(array([[ 1,  2,  1],
        [ 0, -1,  1],
        [ 0,  0,  1]]), array([[ 0.,  0.,  0.],
        [-2.,  0.,  0.],
        [ 1., -1.,  0.]]), array([1, 1, 2]), -1)

In [10]:
U_A, aux_A, c_A, det_A = gauss_elim_nopiv(3, A, b)
print(U_A, aux_A, c_A, det_A)


[[ 1  2  1]
 [ 0 -1  1]
 [ 0  0  1]] [[ 0.  0.  0.]
 [-2.  0.  0.]
 [ 1. -1.  0.]] [1 1 2] -1

In [11]:
x_A = backward_sub(3, U_A, c_A)
print(x_A)


[-3  1  2]

In [12]:
# y_A = forward_sub(3, aux_A, c_A)
# print(y_A)
Função para gerar Matrizes de Hilbert de orden n

2 - Resolver a Matriz de Hilbert com ordem 3, 5 e 10.


In [13]:
def hilbert_gen(n):
    H = np.array([[1.]*n]*n)
    H_b = []
    for b in range(1, n+1): H_b.append(float(b))
    for i in range(0, n):
        for j in range(0, n):
            H[i, j] = (1/((i + 1) + (j + 1) - 1))
    return H, np.array(H_b)

In [14]:
Hil_3, Hil_3_b = hilbert_gen(3)
print(Hil_3, Hil_3_b)
print(np.linalg.solve(Hil_3, Hil_3_b))
print(np.linalg.det(Hil_3))


[[ 1.          0.5         0.33333333]
 [ 0.5         0.33333333  0.25      ]
 [ 0.33333333  0.25        0.2       ]] [ 1.  2.  3.]
[  27. -192.  210.]
0.000462962962963

In [15]:
H3_S, H3_A, H3_B, H3_D = gauss_elim_nopiv(3, Hil_3, Hil_3_b)
backward_sub(3, H3_S, H3_B)


Out[15]:
array([  27., -192.,  210.])

In [16]:
print(gauss_elim_nopiv(3, Hil_3, Hil_3_b))


(array([[ 1.        ,  0.5       ,  0.33333333],
       [ 0.        ,  0.08333333,  0.08333333],
       [ 0.        ,  0.        ,  0.00555556]]), array([[ 0.        ,  0.        ,  0.        ],
       [-0.5       ,  0.        ,  0.        ],
       [-0.33333333, -1.        ,  0.        ]]), array([ 1.        ,  1.5       ,  1.16666667]), 0.00046296296296296005)

In [17]:
Hil_5, Hil_5_b = hilbert_gen(5)
gauss_elim_nopiv(5, Hil_5, Hil_5_b)


Out[17]:
(array([[  1.00000000e+00,   5.00000000e-01,   3.33333333e-01,
           2.50000000e-01,   2.00000000e-01],
        [  0.00000000e+00,   8.33333333e-02,   8.33333333e-02,
           7.50000000e-02,   6.66666667e-02],
        [  0.00000000e+00,   0.00000000e+00,   5.55555556e-03,
           8.33333333e-03,   9.52380952e-03],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           3.57142857e-04,   7.14285714e-04],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           0.00000000e+00,   2.26757370e-05]]),
 array([[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
        [-0.5       ,  0.        ,  0.        ,  0.        ,  0.        ],
        [-0.33333333, -1.        ,  0.        ,  0.        ,  0.        ],
        [-0.25      , -0.9       , -1.5       ,  0.        ,  0.        ],
        [-0.2       , -0.8       , -1.71428571, -2.        ,  0.        ]]),
 array([ 1.        ,  1.5       ,  1.16666667,  0.65      ,  0.3       ]),
 3.7492951325173029e-12)

In [18]:
Hil_10, Hil_10_b = hilbert_gen(10)
print(Hil_10, Hil_10_b)
print(np.linalg.solve(Hil_10, Hil_10_b))
print(np.linalg.det(Hil_10))


[[ 1.          0.5         0.33333333  0.25        0.2         0.16666667
   0.14285714  0.125       0.11111111  0.1       ]
 [ 0.5         0.33333333  0.25        0.2         0.16666667  0.14285714
   0.125       0.11111111  0.1         0.09090909]
 [ 0.33333333  0.25        0.2         0.16666667  0.14285714  0.125
   0.11111111  0.1         0.09090909  0.08333333]
 [ 0.25        0.2         0.16666667  0.14285714  0.125       0.11111111
   0.1         0.09090909  0.08333333  0.07692308]
 [ 0.2         0.16666667  0.14285714  0.125       0.11111111  0.1
   0.09090909  0.08333333  0.07692308  0.07142857]
 [ 0.16666667  0.14285714  0.125       0.11111111  0.1         0.09090909
   0.08333333  0.07692308  0.07142857  0.06666667]
 [ 0.14285714  0.125       0.11111111  0.1         0.09090909  0.08333333
   0.07692308  0.07142857  0.06666667  0.0625    ]
 [ 0.125       0.11111111  0.1         0.09090909  0.08333333  0.07692308
   0.07142857  0.06666667  0.0625      0.05882353]
 [ 0.11111111  0.1         0.09090909  0.08333333  0.07692308  0.07142857
   0.06666667  0.0625      0.05882353  0.05555556]
 [ 0.1         0.09090909  0.08333333  0.07692308  0.07142857  0.06666667
   0.0625      0.05882353  0.05555556  0.05263158]] [  1.   2.   3.   4.   5.   6.   7.   8.   9.  10.]
[ -9.99815326e+02   9.79940514e+04  -2.32814044e+06   2.33001942e+07
  -1.21066246e+08   3.59418665e+08  -6.32245357e+08   6.51054964e+08
  -3.62282609e+08   8.40565858e+07]
2.16440526462e-53

In [19]:
gauss_elim_nopiv(10, Hil_10, Hil_10_b)


Out[19]:
(array([[  1.00000000e+00,   5.00000000e-01,   3.33333333e-01,
           2.50000000e-01,   2.00000000e-01,   1.66666667e-01,
           1.42857143e-01,   1.25000000e-01,   1.11111111e-01,
           1.00000000e-01],
        [  0.00000000e+00,   8.33333333e-02,   8.33333333e-02,
           7.50000000e-02,   6.66666667e-02,   5.95238095e-02,
           5.35714286e-02,   4.86111111e-02,   4.44444444e-02,
           4.09090909e-02],
        [  0.00000000e+00,   0.00000000e+00,   5.55555556e-03,
           8.33333333e-03,   9.52380952e-03,   9.92063492e-03,
           9.92063492e-03,   9.72222222e-03,   9.42760943e-03,
           9.09090909e-03],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           3.57142857e-04,   7.14285714e-04,   9.92063492e-04,
           1.19047619e-03,   1.32575758e-03,   1.41414141e-03,
           1.46853147e-03],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           0.00000000e+00,   2.26757370e-05,   5.66893424e-05,
           9.27643785e-05,   1.26262626e-04,   1.55400155e-04,
           1.79820180e-04],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           0.00000000e+00,   0.00000000e+00,   1.43154905e-06,
           4.29464715e-06,   8.09375809e-06,   1.23333457e-05,
           1.66500166e-05],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           9.00974923e-08,   3.15341224e-07,   6.72727944e-07,
           1.13522841e-06],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           0.00000000e+00,   5.65997065e-09,   2.26398821e-08,
           5.39361899e-08],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           0.00000000e+00,   0.00000000e+00,   3.55136713e-10,
           1.59811304e-09],
        [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
           2.22677000e-11]]),
 array([[  0.        ,   0.        ,   0.        ,   0.        ,
           0.        ,   0.        ,   0.        ,   0.        ,
           0.        ,   0.        ],
        [ -0.5       ,   0.        ,   0.        ,   0.        ,
           0.        ,   0.        ,   0.        ,   0.        ,
           0.        ,   0.        ],
        [ -0.33333333,  -1.        ,   0.        ,   0.        ,
           0.        ,   0.        ,   0.        ,   0.        ,
           0.        ,   0.        ],
        [ -0.25      ,  -0.9       ,  -1.5       ,   0.        ,
           0.        ,   0.        ,   0.        ,   0.        ,
           0.        ,   0.        ],
        [ -0.2       ,  -0.8       ,  -1.71428571,  -2.        ,
           0.        ,   0.        ,   0.        ,   0.        ,
           0.        ,   0.        ],
        [ -0.16666667,  -0.71428571,  -1.78571429,  -2.77777778,
          -2.5       ,   0.        ,   0.        ,   0.        ,
           0.        ,   0.        ],
        [ -0.14285714,  -0.64285714,  -1.78571429,  -3.33333333,
          -4.09090909,  -3.        ,   0.        ,   0.        ,
           0.        ,   0.        ],
        [ -0.125     ,  -0.58333333,  -1.75      ,  -3.71212121,
          -5.56818182,  -5.65384615,  -3.5       ,   0.        ,
           0.        ,   0.        ],
        [ -0.11111111,  -0.53333333,  -1.6969697 ,  -3.95959596,
          -6.85314685,  -8.61538462,  -7.46666668,  -3.99999991,
           0.        ,   0.        ],
        [ -0.1       ,  -0.49090909,  -1.63636364,  -4.11188811,
          -7.93006993, -11.63076923, -12.60000002,  -9.52941157,
          -4.49999391,   0.        ]]),
 array([ 1.        ,  1.5       ,  1.16666667,  0.65      ,  0.3       ,
         0.12301587,  0.0465368 ,  0.01660839,  0.00567211,  0.00187169]),
 2.164480500135833e-53)

O determinante deu 2.164480500135833e-53, que é parecido com o obtido pelo cálculo do determinante utilizando o NumPy

Exercício 3


In [20]:
A = np.array([[1, 2, 1], [2, 3, 3], [-1, -3, 1]])
b = np.array([1, 3, 2])
A, b


Out[20]:
(array([[ 1,  2,  1],
        [ 2,  3,  3],
        [-1, -3,  1]]), array([1, 3, 2]))
3 - Cálculo da matriz inversa utilizando em primeiro lugar a Eliminação de Gauss e Bacward Substitution.

A ideia é resolver n equações para uma matriz quadrada de ordem n, igualando aos vectores.


In [21]:
def matrix_inverse(matrix):
    n = len(matrix)
    inv_matrix = [[] for _ in range(n)]
    for i in range(0, n):
        e = [0, 0, 0]
        e[i] = 1
        U, aux, c, det = gauss_elim_nopiv(n, matrix, e)
        x = backward_sub(n, U, c)
        inv_matrix[i] = x
    inv_matrix = np.array(inv_matrix)
    return inv_matrix.transpose()

In [22]:
A_inv = matrix_inverse(A)
print(A_inv)
print(A.dot(A_inv))


[[-12   5  -3]
 [  5  -2   1]
 [  3  -1   1]]
[[1 0 0]
 [0 1 0]
 [0 0 1]]

Utilizando o NumPy determina-se a inversa de A como:


In [23]:
np.linalg.inv(A)


Out[23]:
array([[-12.,   5.,  -3.],
       [  5.,  -2.,   1.],
       [  3.,  -1.,   1.]])

Exercício 4 - Resolver um sistema para entender a necessidade do pivoting


In [24]:
C = [[-10e-20, 1], [2, 1]]
cs = [1, 0]

In [25]:
solve_partpiv(C, cs)


Out[25]:
array([-3,  2])

In [26]:
x1 = -0.4999975
x2 = 0.999995

In [27]:
C = [[1, 2], [3, 4]]
D = [[1, 1], [3, 2]]